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ABSTRACT 

We discuss the creation of electron-positron cascades in the context of pulsar polar cap acceleration 
models and derive several useful analytic and semi-analytic results for the spatial extent and energy 
response of the cascade. Instead of Monte Carlo simulations, we use an integro-differential equation 
which describes the development of the cascade energy spectrum in one space dimension quite well, 
when it is compared to existing Monte Carlo models. We reduce this full equation to a single integral 
equation, from which we can derive useful results, such as the energy loss between successive generations 
of photons and the spectral index of the response. 

We find that a simple analytic formula represents the pair cascade multiplicity quite well, provided 
that the magnetic field is below 10 12 Gauss, and that an only slightly more complex formula matches 
the numerically-calculated cascade at all other field strengths. 

Using these results, we find that cascades triggered by gamma rays emitted through inverse Compton 
scattering of thermal photons from the neutron star's surface, both resonant and non-resonant, are 
important for the dynamics of the polar cap region in many pulsars. In these objects, the expected 
multiplicity of pairs generated by a single input particle is lower than previously found in cascades 
initiated by curvature emission, frequently being on the order of 10 rather than ~ 1000 as usually quoted. 
Such pulsars also are expected to be less luminous in polar-cap gamma rays than when curvature emission 
triggers the cascade, a topic which will be the subject of a subsequent paper. 

Subject headings: Acceleration of particles — pulsars: general 



1. INTRODUCTION 

The polar-cap pair-production model of pulsar emission 
has remained the chief theory of the region thought to 
give rise to the radiative emission for the past 30 years. 
Particles accelerate from the surface of the neutron star, 
drawn by the induced fields of the rotating magnetosphere. 
As first discussed by Sturrock (1971), these particles emit 
7-rays which pair-produce in the high background mag- 
netic field. These pairs then short out the magnetic field 
(Ruderman & Sutherland 1975), preventing further accel- 
eration. 

The physics of this pair production process has been 
partially explored, with progress ranging from the initial 
generational models of Tademaru (1973), the pair forma- 
tion front structure calculations of Arons & Scharlemann 
(1979); Arons (1981), the detailed Monte Carlo simula- 
tions of Daugherty & Harding (1982), and the more re- 
cent "full-cascade" generational model of Zhang & Hard- 
ing (2000). Our paper describes an improved analyti- 
cal and theoretical understanding of this pair-production 
mechanism, not only reproducing the final 7-ray and pair 
spectra calculated by Daugherty & Harding (1982), but 
detailing the spatial extent of the pair-production process 
itself, giving a more thorough understanding of the pro- 
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cess. 

These straightforward analytic approximations to the 
pair production process were the basis of the model of Hi- 
bschman & Arons (2001), henceforth called Paper I. In 
this paper, we exhibit the full shower theory underlying 
those approximations and find a variety of formulae which 
represent the numerical solutions of the cascade equations, 
thus allowing users of the cascade physics to easily apply 
the results to other problems. We give a full description 
of the pair plasma emergent from pulsar polar caps. We 
reserve discussion of the emergent 7-ray luminosities and 
spectra to a separate paper, as these are of direct obser- 
vational relevance. 



2. GAMMA-RAY OPACITY 

The dominant opacity for pair production over pulsar 
polar caps is the 7-IJ process, wherein a high-energy pho- 
ton interacts with the background magnetic field to form 
an electron-positron pair. The competing process of 7-7 
pair production, wherein a high-energy photon interacts 
with a background X-ray to form a pair, is only com- 
petitive if the background X-ray flux is equivalent to a 
black-body with temperature nearly 10 K, well beyond 
the observed limits on stellar temperatures. 
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According to Erbcr (1966), the opacity for j-B pair- 
production is 

a B (e,ip) = 0.23— - — -sinV'exp (- — I (1) 



Xc B„ 



3x 



where -B 9 is the critical quantum magnetic field, B q = 
e/aX 2 c = 4.41 x 10 13 Gauss, a F = e 2 /?ic « 1/137 is the 
fine structure constant, Xq = h/rac = 3.86 x 10~ n cm is 
the reduced Compton wavelength, B is the local magnetic 
field strength, tp is the pitch angle between the photon 
momentum and the local magnetic field, and 
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where e is the photon energy, in units of mc 2 . For the re- 
mainder of this paper, all energies will be quoted in terms 
of mc 2 , for convenience. 

This expression is accurate provided that esini/> > 1, so 
that the created pair is in a high Landau level, and pro- 
vided that B is small compared to B q ; B < B q /3 suffices. 
For pair production into low Landau levels, the full cross 
section must be used (Daugherty & Harding 1983), and for 
magnetic fields equal to or higher than B qi higher-order 
corrections become important, as discussed in Harding, et 
al. (1997). Most pulsars, however, pair produce into high 
Landau levels for the most significant photons and have 
B < B q /'i. As this work intends to describe the "typical" 
pulsar, we neglect these high-energy and high-field effects. 

2.1. Optical depth 

Because the relativistic primary electron beam follows 
the magnetic field, the primary photons are beamed close 
to parallel to the local field. If the beam electrons have 
a Lorentz factor 7, the primary photons are beamed into 
a cone of opening angle ip ~ 1/7. Typical beam Lorentz 
factors are of order 10 3 or higher, giving an initial pitch 
angle small compared to the pitch angles required for pair- 
production. To a good approximation, then, we can treat 
the photons as if they are injected precisely parallel to the 
field lines. 

As the photons propagate through the magnctosphcrc, 
the pitch angle between the photon momentum and the 
background magnetic field steadily increases. This pitch 
angle is given by sin-0 = |fc x B(r)|, where k is the 
photon momentum direction, and B is the local mag- 
netic field direction. If we assume ip small enough that 
sin^ w ip, the change in ip along the photon path is then 
dip/ds = \kx (fc-V)B(r)|. 

Since ip is small, k remains close to B, so we can substi- 
tute, yielding dip/ds w \B x (B ■ V)JB(r)| = /9(r) -1 , where 
p(r) is the local magnetic field radius of curvature. This 
is accurate to second order in ip. 

For a dipole field, the photon pitch angle depends on 
radius (through second order in the colatitudc, 6) as 



1p = Ipo 
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where r e is the radius of emission and r is the current 
radius, both measured from the center of the dipole, and 
V>oo = r e / Pe , where p e = (4/3)(i?*/^)/ p (r e /i^) 3 / 2 is the 
field line radius of the magnetic colatitude of the field line 
in question at the surface of the star and f p < 1 is a factor 



used to crudely take into account the possibility of non- 
dipolar components to the magnetic field at low altitude. 
Numerically this is Voo = 0Ml{e*/d c )P- l/2 {r/R*) l/2 f P i 
with 9c the colatitudc of the last closed field line, 6c = 
^/R^/Rl and P is the period in seconds. If the dipole 
moment is displaced close to the surface of the star, GR 
bending of the photon paths would change this, but for 
almost all other cases GR perturbations are negligible. 

The magnetic field and radius of curvature may then be 
rewritten as 



B = B P 
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where B e and p e are respectively the magnetic field and 
radius of curvature at the emission point. 

Using these relations to write the optical depth using 
t = ip/ipoo as the integration variable, letting simp ss ip, 
yields 
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dt 
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and e is the photon energy in units of mc 2 . 

The exponential above has a maximum at t — 1/4. 
Physically, this is the point where the increase in opac- 
ity due to the increasing pitch angle balances the decrease 
in opacity due to the steadily declining magnetic field. 

Provided that Z^ is not too small, the integrand is a 
sharply peaked function of t, and can therefore be inte- 
grated using the steepest descents method, yielding 



(8) 

When compared to numerical calculations, this expression 
is accurate to better than 10% in the neighborhood of 
r = 1 , for the range of magnetic fields and radii of curva- 
ture found in pulsars. 

Photons with Too < 1 will escape the magnetosphere, 
while those with r^ > 1 will be absorbed. If Too(e) is 
large, most of the absorption takes place on the leading 
edge of the exponential. Since the opacity is increasing 
exponentially in that regime, we can approximate the in- 
tegral by expanding around the upper endpoint to find 

T(e,t)«A 1 i 3 e- z ~(^ t ) 



and where /(£) has been defined to be t _1 (l — £)~ 3 for 
t < 1/4 and 256/27 for t > 1/4. (This result was previ- 
ously derived in Arons & Scharlemann (1979).) Since the 
opacity decreases rapidly after t = 1/4, the optical depth 
saturates at that point. 
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2.2. Absorption peak 

Given this opacity, we can find where any given pho- 
ton will be absorbed. Due to the exponential dependence 
of the opacity, all photons are effectively absorbed at the 
maximum of a exp(— t), or approximately at r — 1. Using 
(9) for r, we find this peak occurs at 



t a (i-t a y 



(11) 
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where In A = ln(Ai(e)t 3 ). In the limit of small t ai this 
becomes 

*-ia&A- (12) 

Since the opacity saturates past ip a = ^00/4, the minimum 
photon energy required to pair produce is 

32 B i 1 2 n< } \ 

ea = YB^^ mc - (3) 

This energy should be thought of as a critical scaling en- 
ergy, not as the actual minimum energy which will pair 
produce, since the (1 — t) term was neglected. The actual 
minimum energy which will pair produce is 

64 
e mm = -^€ a . (14) 

As long as e > 5e a , the small-t limit (12) is appropriate, 
and we find that a photon will be absorbed after propa- 
gating 

Ar=- — r e . (15) 

4 e y ' 

These expressions can be made tractable by treating In A 

as a constant. Self-consistently evaluating In A at t = 1/8 

then yields (Arons & Scharlemann 1979). 



lnA= 16.2 
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3. SYNCHROTRON EMISSION 



A charged particle with pitch angle i\) with respect to 
the magnetic field emits synchrotron radiation and rapidly 
spirals down to its lowest Landau level. 

The basic rate of synchrotron emission is 

\/3 ac e B 



dN e 



2^- SinV W. K ^ 



(x)dx. (17) 



e/e s 



where e B = B/B q , K 5/3 



is the modified Besscl function of 
order 5/3, and e s is the characteristic synchrotron energy, 
3 



-es7 sin'i/> mc 



(18) 



This corresponds to a total power loss of (Jackson 1975) 

P s = -^eW^ 2 ^mc 2 S -\ (19) 

Due to synchrotron radiation, the Lorentz factor of the 
particle decreases according to 7 = — P s . This energy loss 
occurs over a distance of 1.8 x 10~ 6 7 cm, effectively in- 
stantaneous compared to stellar scales, where we have used 
equation (24) self-consistently to fix jip and let sin-0 s=a ip. 
The final Lorentz factor may be found by noticing that the 
parallel component of the particle's velocity is conserved, 
as can be seen by transforming into a comoving frame. If 
the initial particle has a Lorentz factor of 7^ = (1— /3f )~ 1 / 2 



and is moving at an angle ip with respect to the magnetic 
field, then the final Lorentz factor is 
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(20) 



Integrating (17) over all time gives the number of pho- 
tons generated during the decay. Using (19) to change the 
variable of integration from time to 7 gives 



o AT 

N e (j h ip) = -^—(eBsimpy 1 x 

4e s ( 7 )- 1 / dxK 5/3 (x). (21) 

Changing variables to z = e/e s {^) and re-arranging 
gives 

N e { lu i,) = J^{e BS m4>)- l l 2 e~V 2 [F{ Zl )-F{z f )] 
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where 



F(t) = -J dzz 1 ' 2 J dxK 5/3 (x) 
dxK 5/3 (x)(x^ 2 -t 3 / 2 ) 



(23) 



Equation (22) gives the total number of photons of en- 
ergy e emitted by a particle with initial Lorentz factor 
7i and initial pitch angle ip. Nearly all of these photons 
are emitted in a cone with opening angle ^>, because the 
Lorentz factor of the particle remains high, so the angle 
ip only changes once almost all of the particle energy has 
been radiated away. 

In the case of a particle created by j-B pair produc- 
tion from a photon of initial energy e^, ji = e^/2 and 
^ = V'a(ei)- Defining 
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(24) 



we find Zi = In Ae/e^ and Zf = (1 + a 2 ) In Ae/ej. 

With these substitutions and counting the emission from 
both generated particles, the total number of synchrotron 
photons produced by one incident photon pair producing 

at -0 — ip a is 



N e (e,) 




-3/2 
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We can understand the synchrotron spectrum by exam- 
ining its asymptotic limits. Figure 1 shows the function 
F(t) and its deviation from its zero-point value. The func- 
tion is nearly constant for t < 1, with F(0) — F(t) oc 
t 5 ' e , and exponentially decreases for t < 1. This di- 
vides the synchrotron response into three regimes. If 
e < €i/(l + a 2 ) In A, both z% and Zf are less than one, 
and N e (ei) ex e~ 2 ' 3 ; this is the usual synchrotron spectral 



index. If &;/(! 



1 In A < e < £j/lnA, z t < 1 but Zf > 1, 



and N € (ei) oc e~ 3 ' 2 ; this reflects the decline of the particle 
Lorentz factor. If e > e^/lnA, both Zi and Zf are greater 
than 1, and N e (ei) decays exponentially. 
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Fig. 1. — Synchrotron response function. For t < 1 the function is roughly constant and -F(O) — F(t) is a power law with index 5/6, while 
it exponentially decays at t > 1. This produces the synchrotron spectral indices of —2/3 at low energies and —3/2 at moderate energies, and 
the exponential tail at high energy. 



3.1. Integral equation 

Since particles arc assumed to be beamed along the field 
lines, the specific number intensity at any point is 

7 7 (r, 4>, e) = p(r e )Q 7 (r e ,e)e-^^ (26) 

where r e is a function of r and ip via equation (3), Q 1 is 
the volumetric photon emission rate, p is the radius of cur- 
vature of the field line, and r(r e , r, e) is the optical depth 
attained by a photon of energy e propagating from r e to 
r. The radius of curvature enters since pdtp is the path 
length where the field points in direction dtp. 

The local volumetric rate of pair production events is 
then 

/■OO 

<5 pp (r,e)= / dtpa B (r,ip 7 e)I(r,ip,e). (27) 

Jo 
Using the synchrotron response, equation (25), the syn- 
chrotron emission rate is 

/>00 />00 

Qy,syn{r,e)= de, L N ' e (ej) / di/>aB{r,i>,ei)I(r,ip,ei). 

Jo Jo 

(28) 

Using equation (26) for J 7 , dr = aBpdtp, and the sharp- 
peaked nature of the absorption, we can exactly evaluate 
the ip integral, yielding 



synV, £j 



de l N e (e l ){l - e- T —)Q 7 (f e ,e z ) (29) 



where f e = r e (r,ip), ip is the peak of the ip integral, and 
Tmax is the optical depth attained by a photon of energy 
e propagating from the stellar surface to r. 

Since, by equation (15), a photon emitted at r e trav- 
els at most 0.25r e before pair producing, we can treat the 
cascade as if it occurred on-the-spot (OTS). Photons with 
energies well above e m i n pair-produce in a correspondingly 
shorter distance, so the photons near e m i n set the spatial 
expanse of the cascade. This yields an integral equation 
for the synchrotron spectrum, 

/•OO 

Q- t ,s V n(e)= / <ki(l - e- T ^) X 
Jo 

-K{-)(Q^„ c (€i) + Q 7>8 „„(eO) (30) 



where we have replaced r 7 
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i,, ax with Too and where 
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(31) 



where <f> = (1 + a 2 ). 

This can be easily solved by matrix operator methods. 
If Q ltS yn and Qj. src are represented by vectors Q syn and 
Qsrc, and the integral operator above by a matrix K, the 
solution is simply Q syn = (1 — K)~ 1 KQ src . 

This gives an excellent approximation to the final 7-ray 
spectrum, and a close approximation to the final pair spec- 
trum, as shown in Figure 2. The numerically calculated 
pair spectrum extends to lower energies than the on-the- 
spot (OTS) spectrum, due to the change in magnetic field 
over the course of the cascade. 

3.2. Moment equations 

From the integral equation for the synchrotron photon 
distribution, we find several useful relations connecting the 
moments of the photon spectrum. Using equation (30) to 
calculate successive generations of synchrotron photons, 
we find 



Q { * +1) (e) 



f°° 1 e 
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where Q\ is the i -generation 7-ray spectrum and where 
we have replaced 1 — exp(— t^) with a step function at 

Multiplying by e™ and integrating the above expression 
over all energies gives 



where 
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Q$= dee n Q^(e) 
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i?«= / dee n Q^(e) 



K n = 



dxx n K(x). 



(33) 

(34) 
(35) 
(36) 



R n is the "reduced moment" of the spectrum, neglecting 
energies below e m i„. 

Since K is known, the K n are known. Given a known 

source function, 7?„ may be calculated. The moments of 
successive generations may then be calculated by assuming 



that R n 



Q n so that Qn 



KnQn ■ The limited number 



of generations truncates the sum of successive moments. 
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Fig. 2. — The 7-ray spectrum produced by a single input photon, computed by the OTS model (solid line) and the full numerical model 
(dotted line), evaluated for P = 0.1 s, B = 10 12 Gauss, and e = 10 3 e a = 6.5 X 10 5 mc 2 . 
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Fig. 3. — The pair spectrum produced by a single input photon, computed by the OTS model (solid line) and the full numerical model 
(dotted line), evaluated for P = 0.1 s, B = 10 12 Gauss, and e = 10 3 e a = 6.5 X 10 5 mc 2 . 
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Integrating the kernel, we find that the first two mo 
ments are 

15\/3 



K 



■lnA^ 1 / 2 -!) 



k 1 = i- 4>- 1/2 

In general, the n th moment is 

9V3 2™ In A 1 "" , 
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From ATo and K\ respectively, we can see that the total 
number of photons steadily increases while almost all of 
the energy remains in the photon spectrum, provided that 
4> = (1 + a 2 ) is large. The fraction of energy left in the pair 
spectrum is 1 — K\ = </> -1 ' 2 , which becomes significant for 
values of </> near 1, which occurs at magnetic fields greater 
than 3 x 10 12 Gauss. 

3.3. Generations 

A cascade of pair production then divides into several 
successive generations of photons converting into pairs 
which emit more photons. Following the work of Tade- 
maru (1973) and later papers (Lu, et al. 1994; Wei, et al. 
1997; Zhang & Harding 2000), we find, in our notation, 
that the characteristic synchrotron energies of successive 
generations are related via 



ej+i 



(40) 



In A 

while the particles produced by generation j have final 
Lorentz factors 

7/ , = ^== « £ (41) 



2vT 



or 



where the approximate form assumes a is well greater than 
1, which is true provided that B<3x 10 12 Gauss. These 
correspond exactly to the expressions cited by Zhang & 
Harding (2000), given the equivalence In A — 4/3%. 

From this, it is clear that for each absorbed photon, a 
fraction of the energy, 

fy = 1 - ~J=^ = Ki (42) 

is immediately re-emitted in lower-energy photons, with 
the remainder left in the electron-positron pair. As dis- 
cussed by Zhang & Harding (2000), much of the the en- 
ergy in the pairs is later radiated at energies near 2eb"ff 
via resonant inverse Compton scattering, but this process 
generally does not create further pairs. 

These approximations effectively collapse the entire syn- 
chrotron emission of a created pair into a pulse of photons 
at the peak energy of the synchrotron spectrum. By using 
the moments of the photon spectra, we can derive similar 
approximations, but with wider applicability. 

At each generation, the number of photons increases by 
a factor of Kq, while the total energy in the photons de- 
creases by a factor of K\ , so the average energy of the j th 
generation is 



Q) 



AV-" 



(43) 



where eo is the average energy of the input spectrum. The 
total number of photons in the j generation is 

Nj = N K j Q (44) 

where Nq is the initial number of photons. Provided that 
Zj > e a , most of the photons will convert into pairs. If 
we truncate the generations when the average energy Cj 
reaches e a , we find a power-law relationship between the 
initial energy and the final pair multiplicity, 



M tot (e Q ) ex ( ^ 



(45) 



(46) 



where the power-law exponent is given by 

lnAp 
In A - In AY 

The constant of variation depends on the behavior of the 
lowest-energy photons near e m i n and will be discussed in 
Section 5.2. 

We can apply similar arguments to the spatial extent of 
the cascade. From equation (15), the photons of genera- 
tion j, with mean energy e~j, must propagate a distance 
As = 0.25(c a /ej)R* before pair-producing. Relating the 
multiplicity of generation j, equation (44), to this propa- 
gation distance gives a power-law dependence between the 
multiplicity and the distance, 

M(s) ex s v (47) 

where the power-law exponent, v 1 is identical to that of 
the energy-multiplicity relation, equation (46). 

4. NUMERICAL METHOD 

Since the OTS method cannot spatially resolve the pair 
cascade, we turn to a numerical calculation. Equation (29) 
provides the fundamental structure for the calculation. It 
gives the synchrotron spectrum injected into the magneto- 
sphere at any point, as a function of the photons absorbed 
at that point. 

To calculate the resultant spectrum, we simply evaluate 
the photon spectrum Q(e) at logarithmically spaced points 
in e. The relation between the synchrotron photons out- 
put and the absorbed spectrum at a given radius is then 
a simple matrix multiplication by N e (e'). Due to the vari- 
ation of magnetic field with altitude, this matrix will also 
depend on the altitude. 

We then discretize the radial steps on a variable scale, 
with either purely logarithmic steps, in cases where we are 
primarily concerned with the 7-ray output and the loca- 
tion of the PFF, or a combination of logarithmic steps, up 
to a polar-cap radius, and linear steps beyond, to calculate 
a smooth pair spectrum. Since the final pair energies de- 
pend on the local magnetic field, the pair spectrum is more 
sensitive to the spatial variations in the field, requiring this 
scheme. 

To begin the calculation, at each radial point we in- 
ject an initial spectrum. These spectra are either a delta- 
function source at the lowest radius, for measuring the 
single-photon response, or the expected spectra emitted 
by a single beam particle as it passes through each ra- 
dial bin, for a full physical cascade. In the second case, 
the spectrum is adjusted to prevent binning effects from 
affecting the total power injected. 

Using the expression Too, equation (8), we find which 
photons are absorbed within the magnctosphere, and 



Pair-production multiplicities 



which escape. The escaping photons are simply accumu- 
lated. The absorbed photons are assumed to be absorbed 
at ip(e) — ^> (e), or at a radius of r' = r + p(r)ip a . Each 
photon energy from each radial bin is absorbed at a dif- 
ferent height; these final heights are then simply re-binned 
into the grid. 

For each absorbed photon, equation (25) is used to de- 
termine the synchrotron spectrum injected at the absorp- 
tion point. This process is then iterated, to find the sec- 
ondaries produced by those synchrotron photons, and so 
on, until the energy remaining in synchrotron photons is 
negligible. 

These calculations are streamlined by pre-calculating 
the synchrotron response matrices, N e (ei;r), and the to- 
tal absorption, Too(e;r), at each altitude. In addition, the 
radial absorption matrix, r' '(r; e), is computed for each en- 
ergy, reducing each step of the iteration to a series of ma- 
trix multiplications. 

Figures 4 and 5 show, respectively, the final 7-ray and 
pair spectra for both this model and the Monte Carlo 
model of Daugherty & Harding (1982). The predicted 7- 
ray spectra match extremely well, while the pair spectra 
mismatch at low energies, due to the low resolution of the 
Monte Carlo method in that range, and at high energies, 
due to the approximation in this model of treating In A 
as a constant, rather than allowing it to vary with photon 
energy. 

5. SINGLE-PHOTON RESULTS 

To understand the cascade process itself, we first exam- 
ine the response to a single photon of injected into the 
magnetosphere at the surface of the star. Physically, the 
resultant cascade depends on the energy of the input pho- 
ton, the local magnetic field strength, and the local radius 
of curvature of the field lines. 

5.1. Gamma-ray spectra 

From the discussion of the optical depth, we know that 
only photons emitted at radius r with energies less than 
£mm (r) will escape the magnetosphere. Due to the decline 
of the magnetic field with r, e m i n steadily increases with 
altitude. 

However, the cascade produced by a single injected high- 
energy photon remains localized. If there arc multiple 
generations of pair production, the photon energy must 
be greater than lnAe a , from equation (40), which pair- 
produces in Ar < r e /41nA. Since In A is on the order of 
20, this distance is short compared to r e , and the magnetic 
field remains effectively unchanged. Hence, we may treat 
the cascade process as if it occurred on-the-spot, and the 
energy cut-off remains at e m in(r e ). 

The typical synchrotron response has a power-law ex- 
ponent v = —3/2, extending from energy e^/lnA to 
Ci/(1 + a 2 ) In A, and a power-law exponent of v = —2/3, 
the typical value for synchrotron radiation, for lower ener- 
gies. After processing through multiple generations, this 
translates to an output spectrum with v = —2/3 up to 
£mm/(l + a 2 ) hi A and v = —3/2 from there to e m in- 

Figure 6 shows the gamma-ray spectra produced by vari- 
ous input energies. This demonstrates the expected cut-off 
at e m i n , as well as the transition from a power-law expo- 
nent of —3/2 to —2/3 at e m i„/(l + a 2 ) In A. The shape of 



the spectrum is different at low energies, due to the effects 
of the exponential tails of the opacity and synchrotron 
emission, but it quickly converges to the asymptotic form. 

5.2. Pair multiplicity 

Since the total photon energy at each pair-production 
event is roughly conserved, cascades in low- and mid-field 
pulsars effectively convert high-energy photons into pho- 
tons with energy e m i n . At the penultimate generation of 
an extended cascade, however, a significant fraction of the 
synchrotron photons have e < e m in and escape. From the 
—3/2 synchrotron power law and using the expected mul- 
tiplicity power-law dependence on energy, equation (45), 
we expect a final multiplicity of 

M tot = l 



(48) 



with v as given in equation (46). 

For pulsars with moderate and low fields, B<3x 10 12 , 
the quantity a is large, a negligible fraction of the energy 
of each absorbed photon is left in the generated pair, and 
v w 1 , yielding a simple linear relation between input pho- 
ton energy and final pair multiplicity. 

In the opposite regime, very high fields, a is small, so the 
generated pair retains almost all of the initial photon en- 
ergy, and there is little or no cascade; high-energy photons 
simply transmute into high-7 particles. 

In Figure 7, we plot the numerically calculated energy 
dependence of the total multiplicity along with the theo- 
retical results. The simplest multiplicity formula, equation 
(48) with the exponent v set to 1, works well for magnetic 
fields below 10 12 Gauss, but over-estimates the multiplic- 
ity for higher fields, due to the increase in the amount of 
energy left in the particles in the higher-field cascades. For 
all magnetic fields, the spectrum has three well-defined re- 
gions. If e < e m i n , the initial photon is below threshold, 
and no pairs are produced. If e m i n < e < lnAe m i„, the 
initial photon will pair produce, but all of the secondary 
synchrotron photons have energies below e m i n and create 
no further pairs. Finally, for e > lnAe r „i„, both the initial 
photon and the secondary synchrotron photons may create 
pairs, producing a smoothly increasing multiplicity. 

Figure 8 shows this change in the power-law exponent 
with magnetic field, along with the theoretical result, 
equation (46). The approximation closely matches the nu- 
merical results over the entire dynamic range. 

From equation (47), we expect that the spatial develop- 
ment of the cascade should be a similar power-law in space. 
To find a closer approximation, we use an argument akin 
to that for the final multiplicity. 

Since the first several generations of a high-energy cas- 
cade occur in the high-energy regime where equation (15) 
applies, the cascade develops as if e a were the minimum 
energy, only deviating when the photon energy drops near 
e m i n . Hence, a good approximation to the high-energy 
cascade is 

M(e,s) = l< ' (^( s - s ^y 



27 7^ 
64 4 



Sl<S<Sl . + 64 4" . (49) 

where s is the propagation distance from the emission 
point, si is the point at which the initial photon pair pro- 
duces, si = 0.25(e a /e)r e , v is given in equation (46), and 
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Fig. 4. — Predicted final 7-ray spectrum resulting from injecting a single particle with Lorcntz factor 7 = 6 X 10 7 at the surface of a pulsar 
with a surface magnetic field of 10 12 Gauss and field line radius of curvature of 10 8 cm. The solid line represents the model of this paper, 
while the dotted line represents the results of the Daughcrty & Harding (1982) Monte Carlo model. 
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Fig. 5. — Predicted final pair spectrum resulting from injecting a single particle with Lorentz factor 7 = 6 X 10 7 at the surface of a pulsar 
with a surface magnetic field of 10 12 Gauss and field line radius of curvature of 10 8 cm. The solid line represents the model of this paper, 
while the dotted line represents the results of the Daughcrty & Harding (1982) Monte Carlo model. 
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Fig. 7. — Final pair multiplicity produced by injecting a single photon of varying energy at the stellar surface. The topmost solid line shows 
the theoretical prediction for B = 10 11 Gauss, while the other lines represent the numerical results for different magnetic field strengths: 
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Fig. 8. — Final pair multiplicity power-law exponent, v, where Mtot(t) oc e v . The solid line shows the computed value, while the dotted 
shows the theoretical prediction. 



the factor of 27/64 follows from the truncation of the cas- 
cade at € m in rather than e a . 

Examining the numerical results, shown in Figure 9, we 
found that, indeed, the pair production rate follows this 
rule very well. The total length of the cascade is approxi- 
mately 0.11 i?*, which is the pair production distance ex- 
pected for a photon of the minimum energy capable of 
pair-producing, 64e a /27, according to the tp a model. The 
computed power law, shown in Figure 10 follows the pre- 
dicted values almost exactly. 

5.3. Pair spectra 

Figure 11 shows the variation in the output pair spec- 
trum, for different values of the input energy. The pair 
spectra are characterized by a sharp rise near e a /a, fol- 
lowed by a v — —3/2 power law up to eo/ In A, where eo is 
the energy of the injected photon, followed by an exponen- 
tial decline. The spike corresponding to the input photon 
is also visible. 

Figure 12 shows the variation in the output pair spec- 
trum, for a fixed ratio of input photon energy to e a , 
eoAa = 10 4 , for different values of the magnetic field. 
Here, we notice the similarity of the pair spectra at dif- 
ferent field strengths. The lower limit of the pair spectra 
is at e a /(l + a 2 ) 1 / 2 , which is independent of B, provided 
that a 3> 1. The upper limit is a fixed multiple of this 
value, and so remains constant as well. 

6. SEMI-NUMERICAL MODEL 

The analytic models of Paper I did not take into account 
the variation of the power-law exponent, v, with magnetic 
field. This does not affect the calculation of cither the 
curvature or the resonant ICS cascades, since those mech- 
anisms operate primarily on the primary photons, but the 



decline of v with B limits the effectiveness of non-resonant 
ICS for magnetic fields above 10 12 Gauss. 

Allowing for the change of the power law makes the ana- 
lytic form impractical, but we can quantify the approxima- 
tion by using a simple numerical model. If we still assume 
that all photons emitted by a given emission mechanism 
are emitted at the peak energy and use equation (49 for 
the spatial range of pair-production, we can calculate the 
location of the PFF. 

In practice, we find that this method yields results 
which are nearly indistinguishable from the full numerical 
method, in the regime where the full numerical method 
predicts a finite PFF height. In Figure 13, we show the 
calculated PFF heights from the analytic model and in 
Figure 14 those from the semi-numerical model. 

7. FULL CASCADE 

The same numerical procedure used to calculate the ex- 
pected response to one photon may be applied to a full 
model of a pulsar polar cap. We will ignore any variation 
across the polar cap, concentrating on a typical field line 
instead. 

Given the pair production physics discussed above, a 
polar cap model only requires knowledge of the electric 
potential accelerating particles off the surface of the star, 
the relevant emission mechanisms which inject high-energy 
photons into the magnetosphere, and a model for the tem- 
perature of the star. 

As a model of the accelerating potential, we use the cu- 
bic approximation to the Muslimov & Tsygan (1992) GR 
acceleration derived in Paper I. This model is easily com- 
putable and sufficiently accurate for our purposes. 

The emission mechanisms included are curvature radi- 
ation and inverse Compton scattering, both resonant and 
non- resonant. The spectrum of curvature radiation is well- 
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Fig. 9. — Cumulative pair production as a function of distance, for several values of the magnetic field B and <$> c 



10 14 V. The topmost 



solid line shows the theoretical prediction for B = 10 11 Gauss, while the other lines represent different magnetic field strengths: f X 10 11 , 
3 X 10 11 , 1 X fO 12 , 3 X fO 12 , and f X 10 13 Gauss, respectively. 
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Fig. fO. — Variation of the spatial power-law exponent with B. The solid line is the computed power-law exponent and the dotted line the 
theoretical prediction. 
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Fig. 11. — Final pair spectra produced by injecting a single photon at the stellar surface, evaluated at a B = 10 12 Gauss and <3? ca j> = 10 14 
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Fig. 12. — Final pair spectra produced by injecting a single photon at the stellar surface, evaluated at fixed voltage, <& C ap = 10 14 V, and 
fixed input energy ratio, e = 10 4 e a (B), for various values of the magnetic field. The magnetic fields used arc 1 X 10 , 3 X 10 , 1 X 10 12 , 
3 X 10 12 , and 1 X 10 13 Gauss, from highest peak amplitude to lowest. 
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known, c.f. Jackson (1975), while we use equation (A17) 
for non-resonant scattering, with spatial dilution terms in- 
cluded, and equation (A2) derived from Dermer (1990) for 
resonant scattering. See Appendix A for descriptions of 
the precise emissivities used. 

The stellar temperature is the most uncertain quantity 
in these detailed models. Not only is the precise cool- 
ing rate for neutron stars still a topic of debate, but the 
polar cap is likely heated by beam particles reversed and 
accelerated back down onto the cap. To minimize the un- 
certainties, we simply assume a fixed stellar temperature 
of 10 6 K. 

From the emission mechanisms, we find the Lorentz fac- 
tor of the beam as a function of altitude, and may then 
inject the appropriate emission spectrum at each radial 
bin in our grid. To ensure energy conservation, we slightly 
adjust the peak amplitude of each injected spectrum to 
compensate for the effects of the discrete energy grid. 

We run the model multiple times. First, given an as- 
sumed polar cap temperature, we run it as described above 
to find the pair formation front altitude, namely the alti- 
tude at which n g pairs have been produced per primary. 
At this point the pair plasma is sufficiently dense to short 
out the accelerating electric field. 

We then recompute the Lorentz factor of the beam, with 
acceleration halting at the PFF, inject the new spectra, 
and follow the cascade to produce the final output 7-ray 
and pair spectra. 

A typical output pair spectrum is shown in Figure 15, 
and a typical output 7-ray spectrum in Figure 16. These 
show many of the characteristics of the single-photon re- 
sponse. The pair spectrum shows the rapid rise, v = —3/2 
power-law, and exponential tail, but it also shows a low- 
amplitude tail of high-energy particles, due to the effects 
of inverse Compton scattering. 

For most pulsars, the final pair and 7-ray spectra arc 
similar to those in Figures (15) and (16). 

In Figure 17, we plot the multiplicity as a function of 
pulsar period and cap potential for a fixed stellar tempera- 
ture of 10 6 K. Given that fixed temperature, most pulsars 
have multiplicities ranging from 0.1-100, while most of the 
millisecond pulsars have lower multiplicities. In Figure 18, 
we set the radius of curvature of the field lines equal to 



the stellar radius, in order to model roughly the effects of 
an offset dipolc. Here, most pulsars have multiplicities in 
the range of 10-1000, while the MSPs remain with signif- 
icantly lower multiplicities. However, a parallel offset of 
the dipole moment would also increase the surface field, 
effectively increasing the cap potential. Given the sharp 
gradient in multiplicities at low-P, this could easily bring 
the multiplicities of the MSPs into parity with those of the 
other pulsars. 

8. CONCLUSIONS 

This paper has given several approximations and de- 
scriptions of the pair creation process in pulsars, in the 
hopes that they will be useful for other researchers ex- 
amining these objects. The development of the pair cas- 
cade in space is deceptively simple, a fact which has been 
used in Paper I to obtain several useful results about the 
regimes of dominance of the various emission mechanisms. 
A more detailed model, maintaining the variability of the 
power-law exponent, produces qualitatively identical re- 
sults. Many pulsars have their PFF set by non-resonant 
ICS, with curvature dominating at high potential, and res- 
onant ICS at high magnetic field. 

In addition, due to the effects of NRICS, many pulsars 
seem to operate with comparatively low pair multiplicities. 
The total k is in the range of 1-10 more often than it is in 
the 1000's predicted by curvature models. 

The OTS approximation to the pair cascade describes 
the spectrum of 7-rays and particles produced by the pair 
production process very well, although the OTS pair spec- 
trum underestimates the number of low energy pairs, due 
to the decline of the magnetic field with altitude. 

The full cascade model may be used to predict the 7- 
ray and pair spectra of individual pulsars. Although the 
pair spectrum produced will be modified by RICS, the 
high-energy 7-ray spectrum predicted should be observ- 
able. The 7-ray output from pulsars will be the subject of 
a subsequent paper. 

We wish to acknowledge the assistance of Alice Hard- 
ing, who provided the data for the comparisons of Figures 
4 and 5. 



APPENDIX 

EMISSION MECHANISMS 

For convenience, we describe the forms of the various emission mechanisms used in the numerical calculations. 

Curvature emission 
For curvature emission, we use the standard formula, 

8 2 N C V3a F c 1 [°° 

7uJT~ = T~ "1 2 / dxK 5/3 (x) 

dtde x 3tt A c 7 Je/e c 

where e c = (3/2)Xcpj~ 3 , and K 5 / 3 is the modified Bessel function of order 5/3. 

Resonant ICS 
For resonant ICS, we use the result from Dermer (1990), 

8 2 N r la F c T 7 TA M / e B 

mT, = 2^7^ ^^ l ei ' V 7£B 



(Al) 



(A2) 



where H is the top hat function, defined before. 
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Fig. 15. — Expected pair spectrum for J1952+3252. 
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Fig. 16. — Expected gamma-ray spectrum for J1952+3252, showing the expected v = —3/2 dependence. 
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Final multiplicity 




Fig. 17. — Contour plot of expected multiplicity vs. P, <I>, for a fixed stellar temperature of 10 6 K. 



Final multiplicity (rho = R_*) 




Fig. 18. — Contour plot of expected multiplicity vs. P, <J>, for a fixed stellar temperature of 10 6 K, with a field line radius of curvature, p, 
equal to the stellar radius. 
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Non-resonant ICS 

Although the curvature emission and resonant inverse-Compton scattering spectra are well-represented in the literature, 
the non-resonant ICS spectrum is not, to the accuracy we desire. 

The chief difference between the physical situations considered by Blumcnthal & Gould (1970) and conditions in the 
pulsar magnetosphere is that the background X-rays in pulsars are not isotropic. They originate from the hot surface of 
the star itself and only occupy a range in 9 from 9 = to at most 9 = it. In terms of /i = cos 9, thermal photons are 
emitted with angles from 1 to /i c , where 

^ = -t4=^ ( A3 ) 

V« + z 
where a is the radius of the polar cap, typically a — 9 C R*, and z is the height above the stellar surface. For simplicity, we 
only consider emission along the z-axis, where the incoming photons are symmetric in (f>. The range of /i populated with 
photons is then 

A M =1- Mc = l- / -i— , ■ (A4) 

For the isotropic case, A/i = 2. If A/i < 2, the simplest consequence is that there are fewer photons, by a factor of 
A/i/2, potentially scattering off of the particle. Furthermore, these photons are concentrated near /j, = 1, so the mean 
energy of the photons in the particle rest frame is reduced. 

From (Blumcnthal & Gould 1970), the full non-resonant inverse Compton cross section is 

^■/O - V ( !iV (L + fi _ (i _ „*0\ 8(e' ) (A5) 

dtfti ~2 r °UJ U + e' (1 ^>) d ^ l + (eVmc2)(l+ M ' 1 ) j (A5) 

The primed frame is the particle rest frame, the unprimed frame is fixed to the neutron star itself, and the subscript 1 
indicates the scattered particles and angles. We have made the approximation that all of the incoming photons in the 
particle rest frame are beamed down the z-axis, so that the angle between the incoming and the scattered photon is simply 
equal to the latitude of the scattered photon, ©' = n — 9[. This approximation works very well for large Lorentz factors, 
but is less accurate in the Thomson regime. However, it gives negligible error in the Klcin-Nishina regime. 
The rest-frame scattered number emissivity is then 



jVi.M'i) = 2W|&V dff^f } J'(eV). 



(A6) 



Using the <5-function in the cross section to perform the e' integral yields 

jViX) = \n'rl Q + | - (1 - M?)) 2nJdn'I'(e' ,v,') (A7) 



where 



u l-(e' 1 /mc 2 )(l + M ' 1 ) 
The integral over the source intensity is 



(A8) 



Changing variables to e, using e = 7Cq(1 + P/i'), this becomes 

2n [ dft' I'(e' , M ') = ^\ /'"*" rfe £- 2 /( e , M (e)) (A10) 

The limits of the e integration are set by which energies may be scattered such that they have energy e' in the rest frame. 
The minimum possible energy is e m i n = e(,/7(l — /3/i c ) ~ e o/7^Mj while the maximum is e max — Cq/7(1 — /3) s=s 27£g, 
where the approximations assume that /)»1. 

The source function we adopt is that of a black-body, emitting at temperature T in units of mc 2 , into a range of solid 
angle S/i = 1 — /i c . This corresponds to a specific intensity of 

J(e ' M) = 4^ A | ex P (e/r)-l g ^ Mc ' 1) (A11) 

where if (x, a, 6) = 1 for a < x < b, and otherwise. 
Combining all these components, we find 

^sy^^ - iii/r ^ft (M2) 

1 c £ o r ln ( ±Z *M-Min ) fA 13) 

27^7/3 ^l-exp(- e ^/ 7 rA M ); l JJ 
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The numerator in the log is equivalent to 1 for all energies with significant emission. We then have, 

7v (ei ) = 2.1/ Wl ^^(ei - ^y)i'(eUi) 
1 a\c T f 7 , , , -1 r , , ei 



f v + 4 - (1 - I'?)) '-',, Ml - exp(-4/7Ta„)) 



where we have substituted Tq/)? c = a 2 F /\c- 

If we then assume that ~ 1 so e^/eg = 1 — e i/7: an d evaluate the ei integral, we obtain 
d 2 N c 1 a 2 pC T r _j / x 

'Vi ,, | a „n, (l-eD + ^^-l + Mi 



SiSei 2tt A c 7 4 /?(l - ei) J (1 + /Vi) 2 V 1 - ei 

' ^fi 

^ 2 (l + /?Mi)(l-e"i)rA M 



H ' " exp ( „ 2/1 ; a „ M f 1 1 ,_ VT1A „ ) ) (A14) 



where we have defined e\ = £1/7. 
Changing variables to 



y = -ori , 577^ ,-^a.. . ( A15 ) 



7 2 (l + /3/ii)(l-e"i)rA Ai 
gives, after setting 0=1, 

N(ei) = ^^^^%y / d y (1 + (Mi 2 - 1)(1 - ei) + (1 - ei) 2 ) ln(l - e"«). (A16) 

Approximating this further by replacing fj,' 2 with its average over the range, (n'i), and performing the integral yields 



tf(ei) - ^t^ Jt^TT ^ + «"?> - W 1 - *) + ^ " ^ Li2 ^) 



(A17) 



where 

/ dyln(l-e- ! ')=Li2(e- J/ ) (A18) 



_£l_ 

2 7 2 (l-e- 1 )TA M 



2/min = TT^Ti ,-VnA .. ( Al9 ) 



(l-el)TA M 
where Li2(x) is the dilogarithm function, defined as 



W* = 77 !* (A20) 



Li «(-) = E-^- ( A21 ) 

4 = 1 

The dilogarithm, Li2(x) = x for small x, and Li2(l) = 7r 2 /6. 

The angle average, (/z'] 2 ), clearly lies between and 1. Its value is given by 

„ C dz (z - l) 2 z- 2 (- ln(l - exp(-K/z))) 

^ = A, -2< in < 1 \\\ < A22 > 

J Zo dzz 2 (-ln(l-cxp(-K/z))) 
where 



■■ - , t (A23) 

7 A M T(l- e - 1 ) v ; 

z=l + /Vi ( A24 ) 

z = 1 - w (A25) 

zi = 1 + /3 w 2 (A26) 

Adopting the approximate forms of zq and 21 does not perceptibly change the value of the integral, but reduces (//j 2 ) to 
being a function only of the single parameter k. A good approximation to the value of < /i' 2 > is 

(nf) « 1 - 0.76exp(-(ln K ) 2 /10) (A27) 

For the computation, we simply tabulate (fi'i)- 

The power emitted by the scattering has two asymptotic limits. If "/TA/j, < 1, then emitted power matches the Thomson 
value, 

Pi Th) = |^^7 2 A M 3 T 4 , (A28) 
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while if 7TA/1 > 1, the emitted power is approximately 



p(KN) 



-T 2 A/iln27TA / u 
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(A29) 



12 A c 

which matches the Klein-Nishina limit of Blumenthal & Gould (1970) when A/j, = 2. The transition between the two 
limits is purely a function of 7TA/Z. The emitted power is 



Pn 



1 ^ C T 2 A M /( 7 TA M ) 



2tt A c 
Where / is a numerically calculated function which has asymptotic limits 

f(x) RJ ^X 2 , X«l 

f(x) « ^lo g 2a;, a: » 1 



(A30) 



(A31) 
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